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ABSTRACT. We describe the formation and evolution of spatial and temporal pat- 
terns in cylindrical premixed flames. We consider the cellular regime, Le < 1, where the 
Lewis number Le is the ratio of thermal to mass diffusivity of a deficient component of 
the combustible mixture. A transition from stationary, axisymmetric flames to stationary 
cellular flames is predicted analytically if Le is decreased below a critical value. We present 
the results of numerical computations to show that as Le is further decreased, with all 
other parameters fixed, traveling waves (TWs) along the flame front arise via an infinite- 
period bifurcation which breaks the reflection symmetry of the cellular array. Upon further 
decreasing Le we find the development of different kinds of periodically modulated trav- 
eling waves (MTWs) as well as a branch of quasiperiodically modulated traveling waves 
(QPMTWs). These transitions are accompanied by the development of different spatial 
and temporal symmetries including period doublings and period halvings in appropriate 
fNJ . coordinate systems. We also observe the apparently chaotic temporal behavior of a disor- 

dered cellular pattern involving creation and annihilation of cells. We analytically describe 
the stability of the TW solution near its onset using suitable phase- amplitude equations. 
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Within this framework one of the MTWs can be identified as a localized wave traveling 
O \ through an underlying stationary, spatially periodic structure. 
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1. INTRODUCTION. It has been long known that premixed gaseous combustible 
mixtures do not necessarily burn in a uniform manner. In particular, under certain 
conditions flame fronts can exhibit a cellular structure, characterized by periodic 
arrays of crests along the flame front pointing in the direction of the combustion 
products. The pointed crests are connected by smooth troughs that are convex toward 
the fresh fuel mixture. The temperature is higher at the troughs, which therefore 



appear brighter, and lower at the crests, which therefore appear darker [Q. In many 
instances the transition from laminar to turbulent combustion occurs as a smooth 
flame breaks up into cells, which then undergo transitions leading to increasingly 
complex spatial and temporal behavior. 

In both stationary and propagating cellular flames were observed on a slot 
burner, leading to a linear array of cells, while rotating cellular flames were observed 
on a burner of circular cross section. A variety of cellular patterns have been observed 



for flames stabilized on a circular burner. In [28| various patterns of ordered cellular 
flames were observed. For some of the patterns the cells appeared to be essentially 
stationary, although there was a very small background fluctuation of the intensity, 
the origin of which was not clear. The cells appeared to be reflection symmetric with 
respect to both intensity and shape. In [[31] both rotating and modulated rotating 
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cellular states were observed. The rotating cellular flames were characterized by an 
essentially uniform rotation rate with no visually detectable change in the cellular 
structure as the cells rotated. The modulated rotating state exhibited a modulation 
in cell size, intensity and angular velocity. In both cases the cells appeared to be 
asymmetric in intensity and shape, i.e., the reflection symmetry was broken. 

In |29|] a different type of dynamical state was observed, in which one or more 
members of a cellular array underwent a rapid motion while the other members of 
the array were nearly stationary. In these states, termed hopping states, the hopping 
cells were noticeably asymmetric while the nearly stationary cells appear to be nearly 
symmetric. We note that this type of mode was observed earlier where it was termed 



a square dance mode |25[ . 

In addition to laminar dynamical behavior, chaotic behavior has also been ob- 
served. In [|30] four different manifestations of chaotic cellular flames are character- 
ized. In particular disordered states are observed, characterized by erratic spatial and 
temporal behavior, exhibiting creation and annihilation of cells as well as a chaotic 
variation in intensity, location and cell size. Similar behavior has been observed in 
3§. 

One scenario for the development of cellular flames and the subsequent develop- 
ment of complex spatiotemporal behavior is that of thermo-diffusive cellular insta- 
bilities |37|, |4l], [51] , which are obtained from the diffusional thermal model |42 



in 



which the thermal expansion of the gas is assumed to be weak. The model accounts 
for transport and reaction of heat and one or more components of the combustible 
mixture, as well as advection due to the flow field, which is specified in advance from a 
nonreacting flow calculation. Many of the characteristics of observed cellular flames, 
including complex spatiotemporal flame patterns, can be qualitatively described as 
thermo-diffusive instabilities. Analytical methods have been successful in describing 
primary transitions to cellular flames (e.g. p7j j37j 41 , @)- ^ n some cases higher order 
transitions can be found by analytical methods as well [[4*3] , 13]. In most instances, 
however, secondary and higher order transitions are difficult to find analytically and 
are more readily found by numerically continuing the analytically predicted solution 
branches until transitions are found and then following the new branches. 

In this paper we consider cylindrical flames, in particular adiabatic flames sta- 
bilized by a line source of fuel of strength 2tck. We assume no axial variation and 
restrict attention to an axial cross section. Thus we solve the problem for a fixed axial 
slice and employ polar coordinates within this cross section. We further assume that 
there is a deficient component of the combustible mixture, so that when this compo- 
nent is depleted no further reaction occurs. The model then reduces to two coupled 
reaction diffusion advection equations. It is characteristic of combustion problems 
that the activation energies are large. As a result the reaction terms are important 
only in a thin layer called the reaction zone. In the limit of infinite activation en- 
ergy the reaction zone shrinks to a surface called the flame front across which certain 
jump conditions hold. While we consider the finite activation energy case, for which, 
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strictly speaking, there is no front, we will use this terminology when appropriate. 
Away from the reaction zone the variables change more gradually The cellular pat- 
terns are most pronounced in the narrow reaction zone, where the dynamics are very 
sensitive to the accuracy of the computation. 

The computations build upon and extend analytical and numerical results previ- 
ously obtained for this problem. We first describe the analytical results. There are 
two parameters for this system, k and the Lewis number, Le, the ratio of thermal to 
mass diffusivity. In the infinite activation energy limit this system admits, for all val- 
ues of k and Le, a stationary axisymmetric solution describing a circular flame front 



27\ , f41f . This solution, called the basic solution, is subject to two different classes of 
instabilities depending on whether Le is greater than or less than 1. 

There exists a critical value, Le c \ > 1, such that for Le > Le c \ > 1, the basic solu- 
tion is unstable. The instability arises as two complex conjugate eigenvalues pass into 
the right half plane, thus suggesting that small disturbances evolve to nonstationary 
flames, such as axisymmetric pulsating flames, or flames with traveling or standing 
waves along the flame front. Similar behavior has also been found in other geometries, 
e.g. |37j, |39], [40], [45]]. The regime Le sufficiently greater than 1, describing, e.g., lean, 
heavy hydrocarbon/air mixtures, is often referred to as the pulsating regime. 

Below a second critical value Le C 2 < 1 there exists a k c such that for k > k c the 
basic solution is also unstable. This instability occurs by a pair of real eigenvalues 
crossing into the right half plane, and small perturbations evolve to stationary cellular 



flames [27, 41]. This behavior has also been found in other geometries [35, 37, 43, 51 



The regime Le sufficiently smaller than 1, describing mixtures where the deficient 
component is highly diffusive, for example lean hydrogen/air mixtures or rich, heavy 
hydrocarbon/air mixtures, is referred to as the cellular regime. 



In this paper we consider the cellular regime, Le < 1. The analysis in 127], [|T 



describes the onset of cellular flames as a diffusional thermal instability. However the 
analysis does not describe the subsequent evolution of these cellular flames into more 
complex spatial and temporal patterns and in particular the transition to dynamic 
behavior observed in experiments. These patterns and transitions will be described 
in this paper by numerically continuing the stationary cellular solution branch found 
in |^7J until new branches are found. We then follow these new branches to describe 
more complex spatiotemporal behavior. 

The analysis in |27j ^TJ) identified the roles of the two fundamental parameters 
in promoting diffusional thermal instabilities of cylindrical flames. For the cellular 
regime, Le < 1, the stationary circular flame front can be destabilized by either in- 
creasing k or decreasing Le. In this paper we consider the behavior as Le is decreased 
for a fixed value of k. The behavior of the cellular states as k is increased, with Le 
fixed, is described in ] p~2"f . 

We first describe previous results on cellular solution branches as Le is reduced. 



In |TTJ the behavior of stationary cellular flames was considered. It was shown that 



reducing Le led to a progressive deepening of the cells and more pronounced crests. 
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Only stationary solutions were found for the parameter values used in the computa- 
tion. In H we considered lower values of Le, and found nonstationary 4 cell solutions 
arising via a secondary transition, i.e., a transition from the branch of stationary 4 cell 
flames (S4 branch). These solutions described slowly rotating cellular flames. The ro- 
tation occurred at a constant angular velocity and corresponded to a slowly traveling 
wave (TW) along the flame front. It was shown || that the transition from station- 
ary to rotating behavior occurred as an infinite period, symmetry (parity) breaking 
bifurcation, in which the reflection symmetry of the cells was destroyed. That is, 
while each of the stationary cells was reflection symmetric about its centerline, this 
symmetry was no longer maintained for the TW solutions. We refer to this branch 
as the TW4 branch and will describe it in more detail below. 

The TW4 branch was a pure 4 cell branch, that is, each solution on this branch was 
27r/4 periodic in the angular variable if). Stable mixed mode solution branches were 
also found and described in M. The spatial behavior of these mixed mode solutions 
was that of a 4 cell array in which at any fixed instant of time the cells differed in size 
(e.g., distance between successive minima of the temperature field) and shape. At 
any fixed spatial location the temperature and concentration for these mixed mode 
cellular flames were quasiperiodic in time. Two branches of mixed mode solutions 
were identified. The first branch appeared to be stable over a small window in Le 
near the point where the TW4 branch bifurcates from the S4 branch. The second 
branch was found for values of Le smaller than those for which stable TW4 solutions 
could be found. 

In jnj we presented preliminary results identifying the quasiperiodic mixed mode 
flames as MTWs, in which, after a certain time interval the solution recovers its 
structure along the front except for a phase shift [IE|. Thus the solution becomes 



periodic when viewed in an appropriately chosen rotating coordinate system. In this 
paper we extend the results in ||14|| , identifying additional MTW branches as well as a 
new MTW-like branch in which the cells are quasiperiodically rather than periodically 
modulated. 

We note that in our computations the periodic direction is the polar angle if). 
Thus the period for physically relevant solutions is fixed (2ir) and is not a parameter 
to be varied as in other problems. Consequently the number of cells associated with 
the MTWs, as well as the stability of the various solution branches, is determined by 
parameters related to the gaseous mixture (e.g., k and Le) rather than by an assumed 
period of the solution. For the parameters considered here, the MTW and MTW- 
like branches involve 4 cells rotating about the cylindrical axis. While results are 
presented here for a fixed value of k, we believe that similar results, possibly involving 
a different number of cells, would occur for other values of k. There are complex 
cellular structures and symmetries associated with each of the branches in addition 
to the modulation dynamics. The cellular structures, including symmetries, as well 
as the dynamics will be described in detail here. We also point out the relationship 
of the calculated MTW solutions to experimentally observed flames pi, plj. While 
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our results are for a simplified model of combustion, the striking similarities between 
the computed modes and those observed in experiments indicates that the cellular 
dynamics predicted here does indeed describe those occuring in real flames. 

Specifically, we identify three distinct MTW branches, each differing from the oth- 
ers in the nature of the modulation and the symmetries associated with the modula- 
tion of successive cells. These are referred to as the PMPY (Pushmi-Pullyu) branch, 
the BMTW (breathing MTW) branch and the HPMTW (half-period MTW) branch, 
each of which is described below. In addition, we identify a quasiperiodically modu- 
lated traveling wave (QPMTW) mode, characterized by a quasiperiodic modulation 
involving two apparently independent frequencies. Since there is a third frequency 
associated with the rotation rate, this mode has three apparently incommensurate 
frequencies thus describing a 3-torus. 

We now briefly describe the nature of the solutions on each of the MTW branches. 
Solutions along each branch are composed of 4 cells, which exhibit varied dynamical 
behavior. Solutions along the PMPY branch exhibit a modulation which is strikingly 
different from the other modulations that we observe. Along the PMPY branch the 
solution at any instant of time is composed both of nearly stationary, nearly reflection 
symmetric cells and of nonstationary, asymmetric cells. Motion of the individual cell 
proceeds as a localized wave of asymmetry travels through the underlying stationary, 
symmetric cellular array. This behavior is analogous to the hopping modes observed 



experimentally in |29j and also in |25] where it was termed the "square dance mode" . 
Along the other two MTW branches the rotation rate about the cylindrical axis is 
a periodic (nearly sinusoidal) modulation of the uniform rate characteristic of the 
traveling wave (TW4) branch. The cells also undergo a periodic (nearly sinusoidal) 
expansion and contraction in both size and intensity (maximum temperature within 
the cell). In a frame moving with the uniform rotation rate the cells would appear to 
undergo a breathing motion. The solution is characterized by the symmetry that the 
modulation is the same for each cell, modulo a constant (in time) phase difference be- 
tween any two cells. For one of the branches there is no other apparent symmetry. In 
view of the breathing motion we refer to this branch as the BMTW branch (breathing 
modulated traveling wave). On the other solution branch alternate cells undergo a 
modulation with no phase difference between them, so that the solution is periodic 
over the interval < ip < tt. Since the angular period has been halved, relative 
to the BMTW branch, we refer to this branch as the half-period modulated travel- 
ing wave (HPMTW) branch. Finally, the QPMTW branch differs from the MTW 
branches in that the modulation is quasiperiodic rather than periodic and the motion 
appears to be characterized by three independent frequencies. The solution exhibits 
the symmetry that alternate cells undergo the same modulation with a constant phase 
difference between them. Adjacent cells undergo distinctly different modulations. As 
Le is decreased further the QPMTW branch as well as the HPMTW branch appear 
to lose stability to an attractor which seems to exhibit chaotic temporal behavior 
along with a disordered spatial pattern involving cells undergoing an erratic motion 
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and creation and annihilation of cells in an apparently random fashion. This behavior 
is characteristic of weakly turbulent flames. 

MTWs have been found previously in other parameter regimes in combustion 



as well as in other physical systems. In MTWs were found for nonadiabatic 



cellular flames in the pulsating regime. These modes were shown to be stages in 
the development of chaotic behavior for flames near extinction. MTWs have been 
computed for the Kuramoto-Sivashinsky equation describing combustion and other 



areas of application [18], and have also been observed in non-reacting flows e.g., in 



Taylor- Couette flow experimentally fl32| , and computationally |f2l 



The accuracy of the computed solution is very sensitive to the resolution of the 
reaction zone, the location of which is not known beforehand. As a result, adaptive 
procedures are needed to locate and resolve the reaction zone as the solution evolves 
in time. The problem of obtaining adequate resolution of the reaction zone is further 
complicated by the extremely deep cells which occur as Le is decreased. We employ 
an adaptive pseudo-spectral method which has been used successfully for a variety of 
other problems in combustion (see, e.g., [||, f|, |5|). We describe the model in section 
2 and the numerical method in Section 3. In Section 4 we describe our results in 
detail. In section 5 we provide an analytical description of the transition to the TW4 
branch within the framework of a resonant mode interaction. In addition, we employ 
suitable phase-amplitude equations to discuss its stability as well as the transition to 
the PMPY branch. 



2. MATHEMATICAL MODEL. We consider the problem of a flame stabilized 
by a line source of fuel of strength 2kk. We assume that the reaction is limited by 
a single deficient component and is governed by one step, irreversible Arrhenius ki- 
netics. We denote dimensional quantities by~. The unknowns are the temperature T 
and the concentration C of the deficient component. T u and Tj, are the temperatures 
of the unburned and burned fuel respectively and C u is the unburned value of C . 
Other dimensional quantities are the coefficient of thermal diffusivity A, the activa- 
tion energy E, and the gas constant R. We introduce the nondimensional reduced 
temperature and nondimensional concentration by 

e = (f-f u )/(f b -f u ), c = c/c u . 

The spatial and temporal variables are nondimensionalized by 

_ tU 2 _ £jU 
A A 

where U is the planar adiabatic flame speed for the case of infinite activation energy 
in which the reaction term is replaced by a surface delta function. We employ polar 
coordinates (r,ip), and the nondimensionalized flow velocity due to the fuel source is 
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V = -r where r is the unit radial vector. The equations of the diffusional thermal 



model are [41 



S t = A0-^ + CAexpl^^ ^—), (l: 



r 




AC reCV 

C t = — CA exp 

Le r 

Here A is the Laplacian, a = T u /Tb,N = E/^RTj,), Le is the Lewis number, and 
A = Z 2 / (2Le), where Z = N(l — a) is the Zeldovich number. We note that A, which 
is referred to as the flame speed eigenvalue, depends on the nondimensionalization. 
The value employed above arises from the use of the planar, adiabatic flame velocity 
in the infinite activation energy limit Z — > oo. A different nondimensionalization 
would change the spatial and temporal scales but would not alter the basic patterns 
exhibited by the solution. The boundary conditions are 

C^l, O^Oasr^O, (2) 
C — > 0, 6 — ► 1 as r — > oo. 

In our computations these boundary conditions are imposed at points r\ and r 2 far 
from the reaction zone where combustion occurs. The computed results were found 
to be insensitive to changes in r\ and r%. 

The solution to (|])-(@) has been studied analytically in the limit Z — > oo, and 
Le — 1 = p/Z |27j , In this limit the reaction zone shrinks to a surface r = ty(ijj,t), 



called the flame front. The following stationary, axisymmetric solution exists for all 
values of k and Le: 



e 



f (r//t) K + 0(l/Z), r < re, 
[1, r > re, 

c = i-e + o{i/z), 

* = re, 

and is referred to as the basic solution. The dependence of the solution on p enters 
in the 0(1/Z) terms, not written explicitly here. 

This solution becomes unstable if Le < Le c < 1 and re is sufficiently large, or if 
Le > Le c > 1. In the first case a real double eigenvalue (one corresponding to cos 
the other to sin) crosses from the left half plane into the right half plane and a tran- 
sition from the basic solution to stationary cellular flames occurs [^7], Additional 
transitions have been obtained numerically and are described in Section 4. 



3. NUMERICAL METHOD. We employ an adaptive Chebyshev pseudo-spectral 
method in r that we previously developed p, |5|, together with a Fourier pseudo- 
spectral method in tp. To enhance the resolution of the reaction zone, we adaptively 
transform r in order to minimize a functional monitoring the numerical error. 
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We first briefly discuss pseudo-spectral methods. Detailed descriptions of this 
class of numerical methods can be found in [19|, |33|. We transform the domain for 



the variable r to the interval (—1,1). Here we denote the transformed variable by r 
as well. The solution, u, which generically represents either O or C, is then expanded 
as a sum of basis functions 

j 

u~uj = 53aj-0j-(r). ( 3 ) 

3=0 

In a Fourier method the basis functions <f>j are trigonometric. In a Chebyshev method 
the basis functions <f>j are the Chebyshev polynomials, i.e. <pj{r) = Tj where Tj is the 
j th Chebyshev polynomial, 

Tj(r) = cos(j cos _1 (r)). 

The expansion coefficients aj are obtained from collocation, that is, the function u,j 
is forced to solve the equations at a set of J + 1 collocation points Tj. We employ the 
Gauss-Lobatto points, 

Tj = cos(jn/J), 0<j<J, 

as collocation points. Thus in a pseudo-spectral method the unknowns are the solution 
values at the collocation points. The expansion (Q) is only used to compute spatial 
derivatives. The Fourier method is similar, except that trigonometric functions are 
used as basis functions. 

The major advantage of pseudo-spectral methods over finite difference methods 
is enhanced accuracy for a fixed discretization size. In fact pseudo-spectral meth- 
ods exhibit infinite order accuracy when used to approximate infinitely different iable 
functions, that is, the error can be shown to decrease faster than any inverse power 
of J. These methods are commonly used to approximate problems where a high de- 
gree of accuracy is required, for example in the study of fluid dynamical instabilities, 
(see, e.g., [T9|] ). However, although these methods are highly accurate when used to 



approximate functions which exhibit relatively gradual spatial variation, they have 
difficulties in approximating functions exhibiting localized regions of rapid variation, 
such as the temperature rise and concentration decay across the narrow reaction zone. 
Severe spatial oscillations can occur in approximating rapidly varying functions which 



are not well resolved [|T^, |33|. These oscillations can affect the computed dynamics 
and in certain cases lead to the computation of spurious dynamics (for example in- 
adequately resolved computations can indicate chaotic behavior whereas the exact 
solution is in fact periodic). 

If the location of the reaction zone were known in advance, the oscillations could 
be eliminated by introducing a suitable change of coordinates so that in the new 
coordinate system the solution has a more gradual variation. However, typically the 
location of the reaction zone is not known in advance and is one of the objects of 
the computation. In order to realize the benefits (high accuracy) of the pseudo- 
spectral method in computing solutions which vary rapidly in localized regions of 
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space, we developed an adaptive pseudo-spectral method. This method has proven 
to be effective in computing the rapidly varying solutions which occur in combustion 
problems. The method is described in detail in [[3], H|. The description we give here 
will be brief. 

We introduce a family of coordinate transformations of the form 

r = q(s,a) (4) 

where s is the new computational coordinate and a represents a parameter vector 
which is typically of low dimension. We choose a so that in the new coordinate system 
the solution exhibits a more gradual variation and thus is better approximated by a 
small number of basis functions. Since the behavior of the solution changes during the 
course of the computation, appropriate values of a must be chosen adaptively so as to 
adapt to changes in the solution |3], || . The choice of the coordinate transformations, 
which we discuss below, is motivated by the asymptotic notion of "stretching" a 



layer in the method of matched asymptotic expansions [16, 34] for treating singular 
perturbation problems. 

In order to adaptively choose a coordinate transformation in which the spectral 
approximation is more accurate we employ error measures which are computed for 
each value of a until a minimum is found. In this paper the error measure is the 
functional i 

h{g)= (JjL 2 g) 2 /w(s)d s y , (5) 

where 

wis) = yl — s 2 , L = w(s)—. 

as 

It can be shown that this functional gives an upper bound on the maximum norm of 
the error in approximating a function by its Chebyshev expansion ||. 

A very important feature of the adaptive pseudo-spectral method is the particular 
family of mappings. The unknowns B and C vary gradually except in and near 
the reaction zone where rapid changes occur. We employ a family of mappings, 
introduced in |15|], in which functions with properties similar to these are mapped 



into linear polynomials, which can be approximated using a relatively small value of 
J. 

In order to describe this family of mappings, we let / denote the interval — 1 < 
r < 1 and consider the family of functions 

h(r, ai, a 2 ) = s + tan _1 (o; 1 (r — a 2 ))/X. (6) 

The parameters So and A are determined so that @ maps the interval / univalently 
onto itself, 

so = §q— | , (3 = tan _1 (ai(l + a 2 ))/ tan _1 («i(l - a 2 )), 
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A = tan 1 (ai(l — a 2 ))/(l — Sq). 



When «i is large (the analog of e small in singular perturbation problems), is 
nearly discontinuous with a region of rapid variation occurring near r = a 2 , i.e., there 
is an internal layer at r = a 2 - We observe that other functions, e.g., tanh(z), can be 
used in place of tan _1 (z). 
The inverse of (|6|), 



describes a two parameter family of mappings of / which are one to one and onto, with 
the property that h(q(s, a\, a 2 ), «i, a 2 ) is a linear function. If a,\ and a 2 are properly 
chosen the temperature and concentration profiles will be sufficiently similar to (§) 
that the composite function can be represented by a small number of Chebyshev 
polynomials. The parameter a± is a measure of the rapid rate of change of the 
function, while a 2 is related to the location of the layer, i.e. the region of rapid 
variation. In our method these parameters are determined adaptively by minimizing 
the functional (f|). The equations are integrated in time using a first order splitting 
method described in ||. The timesteps are kept sufficiently small so that there is no 
noticeable effect of temporal integration errors. 

In the computation of cellular flames, the location of the reaction zone depends 
on if). In this case the optimal values of a.\ and a 2 will also depend on if). A two 
dimensional adaptive pseudo-spectral method which allows for this dependence has 
been developed and applied to problems in gasless condensed phase combustion in 
M. The resulting coordinate system is nonseparable, thereby requiring additional 
computation as described in ||. The problem considered here is posed in a coordi- 
nate system in which the Laplacian is separable and it is more efficient to use one 
dimensional coordinate transformations which maintain separability (i.e. ot\ and a 2 
independent of if)), even at the expense of not using values of aci and a 2 which are 
optimal for any particular value of if). We therefore compute the functional (|5]) for 
all angular points and minimize the average, rather than minimizing for each value 
of if). We are nevertheless able to obtain a high degree of resolution of the wrinkled 
reaction zone due to the effectiveness of the family of mappings (0) in concentrating 
resolution in a fixed region. As a result, after the mapping is applied and C are no 
longer rapidly varying in any angular direction although the parameters a± and a 2 
are generally not optimal for any particular direction. 

4. NUMERICAL RESULTS. In this section we present the results of our nu- 
merical simulations of eqs.(|l|) at fixed dimensionless activation energy, N = 20, fixed 
temperature ratio of unburned and burned material, o = 0.615, and fixed flow rate, 
k = 14.8, while the Lewis number Le is varied. The timesteps were typically O(10~ 3 ). 
In all cases 101 Chebyshev collocation points were sufficient to accurately compute 




(7) 
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solutions in the radial direction. We generally used 128 collocation points in the 
angular direction, although in some cases we employed 256. There was virtually 
no effect in increasing the number of collocation points in either independent vari- 
able. In the radial direction the computational domain was taken as r\ < r < T2 
with 7*1 = l,r2 = 41. We found virtually no effect in changing these values. The 
computations presented here were obtained at the NCSA and the NERSC. 

An overview of the different solution branches is given in fig. la where we plot the 
maximum norm of the difference between the computed cellular solution and the (un- 
stable) axisymmetric solution as a function of Le for the different solution branches 
that we have found. We generally concentrate on determining the nature of the tran- 
sitions that occur and on the properties of solutions along different solution branches 
rather than on determining the precise numerical values of the control parameter 
Le at which the transitions occur. We note that we solve the initial value problem, 
marching forward in time until a steady state is achieved, so that in general we only 
compute stable solutions. However unstable solutions can sometimes be computed by 
modifying the computer program so as to impose certain symmetries in the angular 
direction. For example, unstable axisymmetric solutions can be computed by allow- 
ing no variation in if). Unstable 4 mode solutions can be computed by restricting the 
computation to the interval < if) < 27r/4 and imposing periodicity. Unstable reflec- 
tion symmetric solutions can be computed by restricting to the interval < if) < 27r/8 
and imposing reflection symmetry. We have not computed unstable solutions which 
are not stabilized by the imposition of specific symmetries. 

In the regime considered here we find 3 branches of stationary solutions, S4, S5 
and S8, which have, 4, 5 and 8 cells respectively. The cells associated with solutions 
along these branches are reflection symmetric. Only the S4 branch bifurcates stably 
from the basic axisymmetric state; the S5 and S8 branches are unstable at onset. The 
S5 branch becomes stable, however, for smaller Le. This behavior is not unexpected 
since in general only the state closest to the minimum of the neutral curve is stable 
to sideband instabilities at onset, whereas the other states become stable at the 
Eckhaus stability (stability to sideband disturbances) limit. When decreasing Le the 
S4 branch continuously merges with the S8 branch, i.e. the first harmonic becomes 
increasingly stronger whereas the fundamental eventually goes to zero. Thus, the S4 
branch bifurcates from the S8 branch in a pitchfork bifurcation. We will not discuss 
the behavior of solutions that arise when the S5 branch becomes unstable. We also 
did not follow the S8 branch for values of Le below those indicated in fig. la. 

Before reaching the S8 branch, the S4 branch becomes unstable to perturbations 
which break the reflection symmetry (parity-breaking bifurcation) ||. To calculate 
solutions on the S4 branch beyond that point we therefore restrict the solutions to be 
reflection symmetric. At the parity-breaking bifurcation point a branch of traveling 4 
cells, TW4, emanates from the S4 branch. These waves are similar to the rotating cells 
observed in [31] which were also found to be asymmetric while the nonrotating cells 
were symmetric. In the computations the rotation rate and the degree of asymmetry 
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0] go to continuously as the bifurcation point is approached, indicating that this is 
an infinite period bifurcation. Strikingly, the TW4 solutions appear to be unstable at 
onset, although the bifurcation corresponds to a supercritical pitchfork bifurcation. 
Perturbed solutions evolve either to a solution on the S5 branch or, in a small window 
of Le, to a solution on a MTW branch which we refer to as the "Pushmi-Pullyu" 
(PMPY) branch, as discussed below. The relevant instability has the character of a 
sideband instability and can be suppressed by restricting the domain of computation 
to < ip < 27r/4 and imposing periodicity. In the full domain the TW4 solutions 
eventually become stable for smaller Le (Le < 0.37). The instability close to onset 
may explain why the rotation rates of the waves observed in the experiments [[JTJ do 
not approach continuously as the transition point is approached. 

When Le is decreased further, the TW4 solutions become unstable again and 
a branch of modulated waves arises in what appears to be a supercritical Hopf bi- 
furcation. Due to the visual appearance of a breathing of the cells we call this 2 
frequency state a BMTW. For even smaller values of Le two additional, different 
branches of modulated waves arise; the QPMTW, which arises from the BMTW 
branch and appears to be a quasiperiodically modulated wave with three frequencies, 
and the HPMTW branch which bifurcates from the (unstable) TW4 branch. These 
branches are shown in detail in fig. lb. In the following we discuss the properties of 
the modulated waves in detail. 

We first consider the PMPY branch. Fig. 2 shows a contour plot of O at fixed t for 
Le = 0.42. The contours have been chosen to emphasize contours of B corresponding 
to values near 0.98 where the cellular structure of the solution is most pronounced. 
Clearly, none of the 4 cells are exactly alike; two of the cells are nearly reflection sym- 
metric, characteristic of the behavior along the branch S4, while two of the cells are 
not symmetric, characteristic of the behavior along the branch TW4. The dynamics 
of this solution is shown in fig. 3 where 0(r*,ip,t) is shown as a function of i[) for 
increasing values of t. Here and in the following r* is chosen so that the circle r = r* 
is in the reaction zone. As the forward crest (temperature minimum) of a given cell 
(cell 3 in the figure if we number the cells from the left) is jerked clockwise (to the 
left in fig. 3), the cell expands, becomes asymmetric and moves clockwise. This, in 
turn, pulls the cell behind it (cell 4) so that it too expands, becomes asymmetric and 
moves clockwise. At the same time cell 3 pushes the cell ahead of it, compressing it 
and causing it to become nearly symmetric and nearly stationary. This process then 
repeats with cells 1 and 2 etc. in turn and can be viewed as a localized wave of asym- 
metry propagating counter-clockwise through the stationary symmetric cellular array. 



This dynamics suggest that the solution is best described as a "Pushmi-Pullyu" |36 
solution. Its behavior is analogous to that of the hopping modes observed experimen 
tally in |25], |29|. Similar patterns have been observed in directional solidification [^6j 
in directional viscous fingering [47] and in Taylor vortex flow [[52 



Alternatively, it can be seen from fig. 3 that the dynamics of the PMPY solution 
can also be characterized as that of a MTW, since after a discrete time interval r the 



12 



pattern reemerges, though shifted in space, 



G{i/>-cr,t + T) = Q(i(>,t). (8) 

Thus, in a coordinate system rotating uniformly with speed c, ip = ip — ct, the solution 
would be periodic. In addition, the PMPY solutions possess the more restrictive 
symmetry 

6(^ + 27r/4,t + r/4) = 9(^,t), (9) 

which connects adjacent cells. It implies that all 4 cells perform the same dynamics 
although shifted in time. In the context of a discrete or cellular fc-mode system with 
periodicity 2tt, the symmetry 

Q{i> + 2ir/k,t + T/k) = Q$,t), (10) 



referred to as "ponies on a merry-go-round" (POM), is often observed [§ 

The BMTW solutions exhibit the same symmetry ([|) as the PMPY solutions. The 
form of modulation, however, is different. This is shown in fig. 4 which illustrates the 
temporal evolution of Q(r*,ip,t). In order to understand why the TW4 branch does 
not persist stably for lower values of Le, we note that the general effect of decreasing 
Le is to increase the size of the cells, where by the size of the cells we refer to the 
elevated high temperature region between the two successive minima. Decreasing Le 
corresponds to increasing the diffusivity of the reactant relative to the diffusivity of 
heat. As a result the reactant diffuses into a region of higher temperature where a 
higher degree of burning occurs. This in turn generates more heat which then diffuses, 
thus raising the temperature over a larger region, which effectively increases the size of 
the cells. For the values of Le considered, the radial location of the reaction zone does 
not change significantly. Therefore the cells must fit into the circumference 2irr*. As 
Le is decreased, progressively larger cells are forced to fit into the fixed circumference. 
Beyond a certain point a more stable configuration is that illustrated in fig. 4 where 
some cells have contracted while their neighbors have expanded. Analysis of the 
spectrum of the computed solutions indicates that the BMTW branch appears to 
arise due to a sideband instability whereby TW4 solutions are modulated by a mode 
1 modulation. 

The following dynamical behavior is observed. Similar to the PMPY solutions 
the cells undergo an overall rotation with a modulation of cell size and intensity. In 
contrast to the PMPY solutions, the rotation rate never gets close to and all the 
cells are asymmetric at all times. Thus the asymmetry is global rather than localized 
as for the PMPY solutions. In order to illustrate the nature of the modulation, we 
consider the maxima of 0(r*, ip, t). We find that for each fixed t, G attains exactly 
4 maxima for < ip < 2tt which we denote by Bj(i = 1, . . .4). Clearly each such 
maximum is identified with an individual cell and may be thought of as a measure of 
the intensity of the cell. We note that we have accounted for the possibility that any 
particular maximum may well lie between two angular collocation points. We do so 
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by accumulating all of the maxima for which ip is fixed over a time interval and then 
taking the largest of these maxima. In fig. 5a the temporal evolution of 0j is plotted 
for 2 adjacent cells (i.e., i = 1,2). The figure clearly shows both that the modulation 
is periodic and is essentially the same for the two cells except for a constant phase 
shift. Similar properties hold for the other two cells as well. In addition, the phase 
shift is the same between all pairs of adjacent cells. In fig. 5b the angular location 
of each maximum is shown as a function of t. Clearly, in addition to their breathing 
motion the cells undergo a periodic oscillation in position and rotation rate. 

Breathing type MTW solutions, somewhat analogous to the BMTW solutions 
found here, have also been found for nonadiabatic, pulsating (Le > 1) flames near ex- 
tinction fTPf . These MTWs also appeared to arise via a bifurcation from a TW flame, 
specifically a 7 cell TW which was simultaneously destabilized by the subharmonic 



modes 3 and 4. However the MTWs found in [10] exhibited a different symmetry 
from the POM symmetry ([TDD . They exhibited the symmetry referred to as "jumping 
ponies on a merry-go-round" (JPOM) [ [Hfl , 

B$ + l x 2n/k,t + r/k) = 0{$,t). (11) 

Whereas the POM symmetry relates each cell to the cell adjacent to it, this symmetry 
relates cells according to the sequence 1,1 + 1,21 + 1, ...(mod k), effectively skipping 
the I — 1 cells between successive members of the sequence. Here k = 7 corresponds 
to the total number of cells. Thus, for example, if at a certain time cell 1 attains its 
maximum temperature then l/7th of a period later in time cell 1 + I would attain 
its maximum. In an experiment a temperature maximum is observed as a spot of 
maximum brightness. Thus, with I ^ 1 the bright spot would be seen as jumping 
over the I — 1 cells between cells 1 and 1 + /. The POM symmetry is a special case of 
(O) with 1 = 1. We note that the general symmetry (|i~T|) has been discussed in E8 



in the context of rotating fluids. In [|T(| transitions to JPOM MTWs with / = 2, with 
I = 4 and with I = 1 (corresponding to the POM symmetry flTDP) were found. Thus 
MTWs exist in two different parameter regimes within the same model of combustion, 
however the mechanism of instability as well as the symmetries of the resulting MTW, 
is different in the two cases. In both cases there are additional transitions leading to 
chaos, but the nature of the transitions is different for the two regimes. 

When Le is reduced to about Le = 0.21 the BMTW branch becomes unstable and 
a transition to the QPMTW branch occurs, as the symmetry (|9|) is broken. Near the 
transition point the symmetry breaking is most pronounced in the phase differences 
between pairs of adjacent cells. This is shown in fig. 6a where we plot ©(r*, if), t) for 
Le=0.205. It can be clearly seen that cells 1 and 2 (numbered from the left) have 
a small phase difference between them (in fact they appear to be nearly in phase). 
Cells 3 and 4 behave similarly. In contrast, the phase difference between cells 2 and 
3 and between cells 4 and 1 is considerably larger. As Le is further reduced from 
the transition point, the differences in the amplitude of the modulation between the 
odd and even cells becomes more pronounced. This can be seen in fig. 6b where we 
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plot 9(r*,ip,t) for Le=0.180. In addition, each cell undergoes a modulation in size, 
similar to the dynamics of the BMTW solutions. However, adjacent cells now undergo 
different modulations and the 4 cell array splits into two pairs. Within each pair the 
modulations are nearly in phase while there is a large phase difference between the 
two pairs. This behavior has the character of a spatial period doubling. We remark 
that in order to facilitate a comparison between the BMTW, QPMTW and HPMTW 
branches, both the total temporal duration and the time difference between successive 
curves, is the same for figs. 6a and 6b as well as for the analogous figs. 4 and 8. 

While figs. 6a and 6b give the visual appearance of a periodic modulation, upon 
closer examination we find that the amplitude of the modulation is quasiperiodic and 
thus the waves are quasiperiodically modulated. This can be seen in figs. 7 where 
we examine the temporal evolution of the maxima Bj. In figs.7a,b the evolution of 
adjacent and alternate maxima, respectively, is illustrated. For each cell the amplitude 
modulation is quasiperiodic, with a high frequency oscillation modulated by a low 
frequency envelope. Alternate cells have the same envelope, and their modulations 
differ only by a constant phase shift with respect to the envelope. Adjacent cells have 
distinctly different envelopes. In fig. 7c we plot the power spectral density for the 
cell with the larger amplitude modulation (i.e. cell 2 in fig. 6b when ordered from 
the left). In addition to the main peak at uj\ ~ 0.263 we find a low frequency peak 
at cl>2 — 0.016 which corresponds to the slow modulation of the envelope. These two 
frequencies do not appear to be related to each other by a simple rational relation. 
Since this spectrum is calculated from the maximal temperature of one particular 
cell, it represents the dynamics in a frame rotating with that cell. In the laboratory 
frame the rotation rate will add a third frequency to the spectrum. The location of 
the maxima 0, as a function of time is shown in fig. 7d. We note that there is a 
modulation in position (and thus in velocity) and that the modulations are distinctly 
different for adjacent cells. 

Following the QPMTW branch further by decreasing Le, we find a transition to 
solutions exhibiting seemingly random creation and annihilation of cells through phase 
slips. For values of Le near the transition point, the chaotic behavior is characterized 
by intervals where cells are created and annihilated without undergoing any organized 
motion, alternating, in a seemingly random fashion, with intervals in which there 
are 4 cells which undergo a roughly organized motion similar to that of QPMTW 
solutions. This general type of behavior is qualitatively similar to the intermittently 
ordered states observed in [30]. We illustrate this behavior in figs. lOa-b where we 
plot 6(r*, ip, t) for Le=0.165. The figures are sequential and contiguous in t, the 
data is presented in 2 separate figures for clarity. Note that although a transition to 
chaos can generally be expected from a branch of quasiperiodic solutions with three 
independent frequencies [fTT , here the transition may not be via the quasiperiodic 



route as a result of symmetries. This is also supported by the fact that phase slips 
are very important in the chaotic dynamics so that the qualitative features of the 
chaotic attractor strongly differ from those of the quasiperiodic attractor. 
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We also find a narrow window in Le in which a fourth branch of stable modulated 
waves exists. The behavior of solutions along this branch is illustrated in fig. 8 where 
Q(r*,i[),t) is plotted for Le = 0.175. Again, the shape of the cells is modulated in 
time. In contrast to the BMTW and QPMTW solutions, however, alternate cells are 
in phase, i.e. the solution has the discrete translation symmetry 

e(^ + 27r/2,f) = e(^*), 

which does not involve any phase shifts. Thus the spatial period is only half the 
system size, as compared to the BMTW solutions which only have the period 2ir. We 
therefore call this new branch a half-period MTW (HPMTW). 

Our results indicate that the HPMTW branch bifurcates from the TW4 branch. 
Since the TW4 solutions have the spatial symmetry 

e(^ + 2vr/4,t) = 0(^,t), (12) 

the bifurcation corresponds to a pure period doubling in space. For comparison, 
the analogous transition from the TW4 branch to the BMTW branch constitutes a 
period quadrupling. The modulation of the temperature maxima is strictly periodic 
as can be seen from fig. 9a where the temperature maxima for two adjacent cells are 
plotted. The angular locations of the 4 maxima are plotted in fig. 9b as a function of 
time. Note that the rotation direction is reversed from the previous cases. This is an 
effect of the initial conditions. For each rotating solution, a solution rotating in the 
opposite direction can easily be obtained by simply reflecting the angular coordinate. 
We also observe from fig. 9b that the oscillation in peak location is exactly in phase 
for alternate cells. 

The HPMTW branch can be obtained by following the (unstable) TW4 branch, 
imposing periodicity over the interval < ip < it. At onset the HPMTW branch 
is unstable when considered over the full angular interval < ip < 2n. This is 
expected since it bifurcates from the TW4 branch which is itself unstable for this 
range of Le. From the computation in the restricted domain we have observed that 
the subharmonic mode (mode 2) appears to grow continuously from as Le is reduced 
thus indicating that the HPMTW branch bifurcates from the TW4 branch. In the 
full domain the HPMTW branch becomes stable for lower values of Le. Both in the 
restricted and the full domain the HPMTW branch loses stability to an apparently 
chaotic attractor for Le approximately 0.165. Thus, its range of stability extends 
beyond that of the QPMTW branch and at Le = 0.170 the HPMTW branch coexists 
stably with an apparently chaotic solution arising from an instability of the QPMTW 
solutions for Le sufficiently small. 

We note that the two successive bifurcations of modulated waves from the TW4 
branch (to the BMTW and HPMTW branches, respectively) suggests the possibility 
of a third bifurcation, to a branch which, extending our notation might be called 
a quarter period MTW, in which each cell is modulated but all cells are exactly in 
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phase and satisfy the symmetry (|T2|). In fact the modulated wave reported in |3T| 
appears to satisfy an analogous symmetry (the number of cells was different). An 
investigation of the TW4 branch for smaller values of Le has failed to find such a 
branch although it is quite possible that such a branch exists for other parameter 
values. 



5. ANALYSIS. In this section we describe two approaches which yield insight into 
the connection between the S4, S8 and TW4 branches as well as between the TW4 
and PMPY branches, respectively. The merging of the S4 branch with the S8 branch 
suggests considering the resonant interaction of a mode with amplitude A and wave 
number q and a mode with amplitude B and wavenumber 2q corresponding to S4 and 
S8, respectively. In the vicinity of a codimension-2 point in parameter space at which 
both these modes simultaneously destabilize the axisymmetric state the dynamics of 
the system restricted to 27r/4 periodic solutions can be described by two complex 
amplitude equations 

d T A = inA + Cl A*B + ai A\A\ 2 + b 1 A\B\ 2 , (13) 
d T B = u 2 B + c 2 A 2 + b 2 B\A\ 2 + a 2 B\B\ 2 , 

where * denotes complex conjugate, and T is a slow time variable. The coefficients a^, 
bi, and c, are 0(1), and n t and /i 2 are the control parameters. These equations, which 
have been studied extensively, e.g. |l], |23], 46 , exhibit a parity-breaking transition 



from a stationary q mode branch to a TW branch before the q mode branch merges 
with the 2q mode branch, just as in the numerical computations of cellular flames 
considered here. Though these simulations are not close to such a codimension 2 point, 
calculations in directional solidification |^9[ and in Taylor vortex flow pO] showed 



that the regime of existence of the TW branch can extend far beyond this special 
point in parameter space. Thus the amplitude equations ( |T3"D for the q : 2q mode 
interaction appear to describe transitions from S4 to TW4 and to S8, respectively. 
The amplitude equations also describe MTW solutions, but these are different from 
the MTW solutions described in this paper. This is due to the fact that eqs.(|T3"D are 
restricted to a description of the behavior of four identical cells, while in our MTW 
solutions PMPY, BMTW and HPMTW as well as in our QPMTW solution, the four 
cells do not behave identically. 

To understand the stability behavior of solutions on the TW4 branch as well as 
the relation of the TW4 branch to the PMPY branch, it is useful to take a different 
approach. In the vicinity of a parity-breaking instability the dynamics can be de- 
scribed by two coupled real equations for the amplitude A of the small antisymmetric 
part of the solution, which measures the extent to which the reflection symmetry is 
broken, and the local phase <fi of the solution, which gives the instantaneous position 
of each cell |22| . In the context of the present computation the symmetric solution 
S4 corresponds to A = and the asymmetric TW4-solution to A = Aq = const. In 
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order to describe the PMPY-solution one would have to allow the asymmetric ampli- 
tude to depend on the spatial coordinate. This is systematically possible only in the 
limit that the spatial variations occur on a long length scale. We therefore imagine 
introducing the radius R of the circular flame front as an additional parameter which 
can be changed by changing the flow rate k. In the limit R — > oo the introduction of 
a slow angular coordinate X = eip and a slow time variable T = et, e -C 1, on which 
A and (ft are allowed to depend, is justified and one can find evolution equations for A 
and (ft. One may hope, that these equations then also describe the system with finite 
R in a qualitative way. Strikingly, this is indeed the case as is shown below. 

To derive the evolution equations for the amplitude A of the asymmetric part and 
the local phase (ft we expand the physical variable, say u, as 

u = u s {(ft) + eA(X, T)u A {(ft) + h.o.t., 



and use symmetry and scaling arguments to obtain the evolution equations [20, 22 



d T (ft = A, (14) 
d T A = X 1 d X (ftA + b 2 d X x<ft 

+e {(A + X 2 (d X (ft) 2 )A + d 2 d xx A - A 3 + gAd x A + b 21 d X (ftd XX (ft} . 

Note that the equations explicitly contain e which signifies that, strictly speaking, two 
separate slow time scales may be called for. Instead, we consider the reconstituted 
equations and keep track of the relative order of the terms as signified by e. Extended 
traveling wave solutions, which we associate with TW4 solutions, are given by 

A 2 = X = Al (ft = uT (15) 

with uj = A . In general, (ft could also have a linear dependence on X, but it can be 
absorbed into the coefficients Ai, A2 and b 2 i- To determine the linear stability of this 
solution we write it as 

A = A + 5A 1 e ipx+at , (ft = (fto + S(ft 1 e ipx+at . (16) 

Inserting into ( Ti~4[) and linearizing in 5 1 yields the dispersion relation 

a 2 + ea (2A 2 + d 2 p 2 - ipgA ^j + b 2 p 2 - ipAoXt = 0. (17) 

Using the Hurwitz criterion, the TW solutions are stable iff 

d 2 > (18) 

and 

H{p, A*) = b 2 d 2 2 p 4 + p 2 d 2 A 2 (4b 2 + \ ig ) 

+ 2X l9 A 4 + Ab 2 A A - XjA 2 e- 2 > 0. (19) 



18 



Thus, stability with respect to perturbations with large wave numbers p (short waves) 
requires 62 > 0, which is also a condition for the stationary structure described by 
A = to be Eckhaus-stable. If 62 < 0, a fourth derivative of <fi must be included 
in fll4]), as in ( p5|) below, in order to obtain a well posed problem. Strikingly, if e is 
small, traveling wave solutions are unstable to long wavelength perturbations due to 
the term involving Ai [ 20| . That is, for e and p small, (pi]) cannot be satisfied. Thus, 



in an infinite system, which allows p — > 0, i.e., allows very long- wave perturbations, 
the traveling waves are in general unstable at onset if Ai 7^ 0. Since Ai determines the 
slope of the neutral stability curve corresponding to the secondary bifurcation, this 
stability behavior is similar to that of stationary structures which arise in a primary 
bifurcation; there, only the structure with the critical wavenumber (at the minimum 
of the neutral curve, where the analog of A 1 is 0) is stable at onset with respect to 
the Eckhaus instability. States with other wave numbers can become stable further 
above threshold. By analogy, we might expect that the TW branch could also gain 
stability at larger amplitudes, even if Ai 7^ 0. This can be seen by assuming Ai to be 
small, 

Ai = eAi. (20) 

Since H is strictly increasing with p 2 , the TW solutions become stable in an infinite 
system when H{p = 0) = 0, i.e. for 

A2 ° 2 < 21 > 

The TW branch A 2 = Ao arises through a supercritical pitchfork bifurcation. In a 
finite system with periodic boundary conditions the bifurcation is also a pitchfork and 
we therefore expect the TW solutions to be stable even at onset, since only a single real 
eigenvalue passes through 0. The eigenvalue corresponding to the translation mode 
is identically zero and, since the spectrum is discrete, the remaining eigenvalues are 
in the left half plane, bounded away from the imaginary axis. This result can be 
recovered from ( |T9"D by limiting p to p > p c = 2tt/L where L is the length of the 
system. Then for small e the TW branch is stable in two windows given by 

(22) 
(23) 

The TW branch becomes stable for all amplitudes, i.e., the window of instability 
disappears, if 

V < ™g£ (24) 
Based on (20) we might therefore expect that the TW solution is stable at onset. 
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In contrast to this argument, in the numerical computation of ([I]) the TW4 branch 
seems to be unstable as close to onset as could be investigated. In addition, the sta- 
tionary S4 branch appears to be stable all the way to the parity-breaking bifurcation, 
beyond which, presumably due to the Eckhaus instability, perturbations evolve to 
the S5 branch. We therefore consider the case that the parity-breaking instability 
and the Eckhaus instability (b 2 = 0) occur extremely close to one another. Such a 
situation occurred in Taylor vortex flow |50"| . Thus we extend equations (0) to the 
case b 2 — e 2 /^ and allow f3 to have either sign depending on which bifurcation occurs 
first. This leads to 

d T <p = A, (25) 
d T A = X^x^A 

+e f (A + A 2 {d x <pf)A + d 2 d 2 x A - A 3 + gAd x A + b 21 d x <J)d x c 



+e 2 {(3d 2 x <p - 6 4 0* « 
Note that we assume 64 > for well posedness. Taking Ai = e 2 A 1 leads to 

H = (4((3 + b 4 p 2 ) + 2A ig )A 4 

+ (d 2 p 2 (A(3 + Ab 4 p 2 + A l9 ) - A 2 ) A 2 + d 2 2 p 4 (p + hp 2 ) . (26) 

Thus, for small amplitudes the TW solutions are unstable to long-wavelength pertur- 
bations p 2 < —6 4 //3, as in the stationary structure. However, for larger amplitudes 
the unstable TW solutions become stable if 

K x g > -2(3 > 0. (27) 

Thus, with (3 < the above scenario describes a TW branch which is unstable at onset 
but stabilizes at some distance from onset, as found in the numerical computation of 

The coupled phase-amplitude equations (|25| ) also give insight into the PMPY 



solution discussed in section 4. The PMPY solutions can be described as a localized 
wave of asymmetry traveling through the underlying symmetric array of stationary 
cells. For the case b 2 > 0, for which the stationary cells are Eckhaus-stable, it has 
been shown that (|T4D supports stable solutions of that kind [20 , [50] • Their existence 



can be seen easily in the special case g = b 2 \ = \ 2 = |50| . After introducing q = d x 



and going into a frame moving with velocity v, (|T4D can be solved for q, as 

q = -A/v + goo. (28) 

For stationary solutions this leads to an 'equation of motion for a particle with mass 
d 2 \ 

d 2 d\A +(v-^\ d x A = -d A U{A) (29) 
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in the 'potential' 



U(A) = ^(eA + \ iqoo )A 2 - ^A 3 - '-A'. (30) 



Localized solutions correspond to homoclinic orbits of ( P§D which start at the max- 
imum of the potential at A = for X = — oo, reach a turning point at A max and 
return to A = for X = oo. Such orbits exist if the 'friction' v — — vanishes and 

V 

if A > A c = — Ai/e (2Ai/(9e6 2 ) + q^). The latter condition ensures that the value 
of the potential at its second maximum, which occurs at nonzero A, is positive and 
therefore above the value at the turning point A max , U(A max ) = 0. 

Since the numerical calculations in section 4 strongly suggest that the stationary 
structure becomes Eckhaus-unstable very close to the parity-breaking bifurcation, the 
analysis of (|i~4]) is not sufficient. Applying a similar analysis to the extended eqs.([25|) 
leads to an integral condition for the determination of v, given by 

6621 , \ 4 a A i { & 2 e&21<7ocA o , . e 2 &4 ~ 3 ,\ , v 

+ eg \ Ad x A + [v d x A H d x A dX 

J \ v v ) v J 

= 0, (31) 

which replaces the condition v — b^jv = 0. This condition expresses the statement 
that there is no change in the total 'energy' along the homoclinic orbit. Instead of 
treating this condition together with the equation of motion by perturbation methods 
pi| , we solve fl2"5p numerically, which, in addition to showing the existence of such 
localized solutions, will also show their stability, at least for the given parameter 
values. 

Localized waves behave like a symmetric stationary state (A = 0) for X — > ±oo, 
i.e. for q — > q^, and like an asymmetric traveling wave (A 2 = X Q ) in a localized region. 
From (2"Bp with A 2 = X Q the wave number in the asymmetric region is q as ^ q^. We 




would like to determine parameter values for use in the numerical solution of (|25|) , 
to exhibit stable localized waves. Thus we would like to determine conditions on the 
parameter values to insure that the symmetric stationary state is stable for X — > ±oo 
(<? — > <?oo), and that the asymmetric traveling wave is stable for X in the asymmetric 
region (q = q as ). It is simpler to determine conditions that insure that the symmetric 
stationary and the asymmetric traveling wave solutions, each defined on the entire 
X interval, are simultaneously stable. These conditions then provide an estimate for 
the parameter values to be used in (^5|) . For the stability of the symmetric stationary 
state we require 

b 2 + b 2iqoo > 0, d 2 > 0, A 1(?00 + e(A + A 2 g^) < 0, (32) 

where the first condition corresponds to the condition 6 2 > given above, with the 
role of 6 2 played by 6 2 + fr^oo, and the third condition follows from the requirement 
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that the growth rate (coefficient of the term linear in A in (|25|)) is negative. For the 
stability of the asymmetric traveling wave we choose g to satisfy 



A l9 > -2(6 a + b 21 q as ) > 0, (33) 

which is the analog of (0), where now b 2 + b 2 iq as plays the role of b 2 . Thus, though 
the stationary state at the local wave number q as would be Eckhaus-unstable since 
b 2 + b 2 \q as < 0, the traveling wave can be stable due to the presence of the nonlinear 
gradient term gAdxA. Fig. 11a shows the temporal evolution of A for such a stable 
localized wave with A = —0.1, Ai = —3, A 2 = 0, d 2 = 1, b 2 = 0, 6 2 i = 0.1, 6 4 = 1 
and g = —0.5. In fig. lib the amplitude and the wave number are shown for a 
representative time. 



Unfortunately, the coefficients for ([14]) are not known for the system of interest 



here. Thus, a quantitative comparison is not possible. The above calculations show, 
however, that these equations may very well describe the behavior of the cellular 
flames in the vicinity of the parity-breaking bifurcation. 
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Figure Captions 

Figure la. Solution branches as a function of Le. 

Figure lb. Solution branches over restricted region in Le. 

Figure 2. Contour plot for 6 for solution on PMPY branch, Le = 0.42. 

Figure 3. G(r*,ip,t) for solution on PMPY branch, Le = 0.42, r* = 13.5. Time t 
(angle ip) increases along vertical (horizontal) axis. 

Figure 4. Q(r*,t/j,t) for solution on BMTW branch, Le = 0.23, r* = 10.4. Time t 
(angle ip) increases along vertical (horizontal) axis. 

Figure 5a. Temperature of two adjacent maxima (peaks) vs. t, for solution on BMTW 
branch, Le = 0.27, r* = 11.158. 

Figure 5b. Angular location of temperature maxima (peak) vs. t, for solution on 
BMTW branch, Le = 0.27, r* = 11.158. 

Figure 6a. Q(n,ijj,t) for QPMTW solution , Le = 0.205, r* = 10.2. Time t (angle 
ip) increases along vertical (horizontal) axis. 

Figure 6b. Q(r*,ijj,t) for QPMTW solution , Le = 0.185, r* = 10.0. Time t (angle 
ip) increases along vertical (horizontal) axis. 

Figure 7a. Temperature of two adjacent maxima (peaks) vs. t, for solution on 
QPMTW branch, Le = 0.185, r* = 9.842. 

Figure 7b. Temperature of two alternate maxima (peaks) vs. t, for solution on 
QPMTW branch, Le = 0.185, r* = 9.842. 

Figure 7c. PSD corresponding to Fig. 7a. 

Figure 7d. Angular location of temperature maxima (peaks) vs. t, for solution on 
QPMTW branch, Le = 0.185, r* = 9.842. 

Figure 8. 9(r„^,t) for solution on HPMTW branch, Le = 0.175, r* = 9.8. Time t 
(angle ip) increases along vertical (horizontal) axis. 

Figure 9a. Temperature of two adjacent maxima (peaks) vs. t, for solution on 
HPMTW branch, Le = 0.175, r* = 9.346. 

Figure 9b. Angular location of temperature maxima (peak) vs. t, for solution on 
HPMTW branch, Le = 0.175, r* = 9.346. 

Figure 10a. 0(r*, ip, t) for solution apparently exhibiting apparently chaotic temporal 
behavior with a disordered spatial structure , Le = 0.165, r* = 9.6. Time t (angle ip) 
increases along vertical (horizontal) axis. 
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Figure 10b. 0(r*, ip, t) for solution apparently exhibiting apparently chaotic temporal 
behavior with a disordered spatial structure , Le = 0.165, r* = 9.6 Time t (angle ip) 
increases along vertical (horizontal) axis. Continuation of fig. 10a. 

Figure 11. 

a) Temporal evolution of asymmetry amplitude A according to phase-amplitude equa- 
tions fl2~5]) for A = —0.1, Ai = —3, A 2 = 0, d 2 = 1, b 2 = 0, b 2 i = 0.1, b A = 1 and 
<? = -0.5. 

b) Representative stable localized wave at same parameters as a). 
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